MODULE OUTVAR
!
! Verwendungszweck:
!
!      Enthaelt alle Ausgaberoutinen sowie Auswerteroutinen zum berechnen wichtiger Groessen
!
!          OUTPUT        : Main Routine - contains the following subroutines
!            CALC_ETOT   : Berechnung von Gesamt-drehimpuls, -masse und -energie
!            WRITE_TPD   : Schreiben des .tpd Files;  R, u, L, T  von der Photosphaere und vom Aussenrand
!            WRITE_DAT   : Schreiben des .dat Files;  komplettes Modell zu einem Zeitpunkt
!
!
!    BEMERKUNGEN
!
!     - alle Routinen in diesem module beziehen sich auf X und die
!       dazugehoerigen Hilfsvariablen und Materialfunktionen
!
!     - waehrend dem output sind jedoch XA und X eh fast immer ident: entweder nach einer konvergenz in 'newton'
!       oder nach 'getready' (einzige Ausnahme: nach einem divergenten Abbruch in 'time').
!       Insofern ist es fast egal ob wir hier X oder XA verwenden...
!
!     - Die Hilfsvariablen beziehen sich i.a. auf den Variablensatz X (bis auf einige 'alte' Hilfsvariable mit XA).
!       Insofern muss man konsistenter weise generell den "neuen" Variablensatz X fuer den Output heranziehen.
!
!     - wir brauchen hier an viellerlei stellen einen aktuellen satz hilfsvariable:
!       Das ist kein kein problem, da vor dem Aufruf von 'output' entweder 'newton' gelaufen ist
!       oder - beim beginn der Rechnung - 'getready';
!       Lediglich im Falle eines divergenten Abbruchs  in 'time' (d.h. nach dem Kriterium tst < tstmin)
!       koennen die Hilfsvariablen und Materialfunktionen um einen Iterationsschritt veraltet sein.
!
!
!-----------------------------------------------------------------------
! Interne Variable:
!
!     L_tot     : Gesamtdrehimpuls des aktuellen Modells
!     L_in_sum  : Integral des Drehimpulstransports ueber den Innenrand in allen bisherigen Zeitschritten
!     L_sum : Integral des Drehimpulstransports ueber den Aussenrand in allen bisherigen Zeitschritten
!     L_err     : relativer Fehler in der Drehimpulserhaltung waehrend eines Zeitschrittes
!
!     m_tot     : Gesamtmasse des aktuellen Modells
!     m_in_sum  : Integral des Massentransports ueber den Innenrand in allen bisherigen Zeitschritten
!     m_out_sum : Integral des Massentransports ueber den Aussenrand in allen bisherigen Zeitschritten
!
!     accr_in   : Akkretionsrate ueber den Innenrand
!     accr_out  : Akkretionsrate ueber den Innenrand
!
!     e_tot     : Gesamtenergie des aktuellen Modells
!     twistmax  : maximale Verdrehung durch die Rotation der Scheibe waehrend eines Zeitschrittes
!
!     lum       : Gesamtleuchtkraft der Scheibe
!     T_eff     : Leuchtkraft-gemittelte Effektivtemperatur der Scheibe
!
!     CFLnumber : Verhaeltnis von Advektionsgeschw. zu Zellgroesse
!
!     save_dat  : Jedes "save_dat"-te Modell wird als .dat-file gespeichert
!     save_tpd  : Nach jedem "save_tpd"-ten Zeitschritt erfolgt eine Ausgabe ins .tpd-file
!
!-----------------------------------------------------------------------

      use config, only : rkind, np

      implicit none

      private
      public :: OUTPUT, save_tpd, save_dat

      real(kind=rkind)    :: L_tot, L_in_sum, L_out_sum, L_err
      real(kind=rkind)    :: m_tot, m_in_sum, m_out_sum
      real(kind=rkind)    :: accr_in, accr_out, e_tot, twistmax
      real(kind=rkind)    :: lum, T_eff

      real(kind=rkind)    :: CFLnumber = 0.

      integer             :: save_tpd, save_dat

CONTAINS




SUBROUTINE OUTPUT(status)
!
! Verwendungszweck:
!
!      Aufrufen der Auswerte & Ausgabe-programme
!
!        'status' kann drei werte annehmen:   ini ... beim ersten aufruf
!                                             run ... waehrend der rechnung
!                                             fin ... beim letzten aufruf nach beendigung der rechnung
!
!-----------------------------------------------------------------------
! Interne Variable:
!
!     k, l, m   : Zaehlvariable bis zum naechsten Output
!
!-----------------------------------------------------------------------

      use global, only : zz, loglevel
      use physco, only : z0

      implicit none

      integer                      :: k=0, l=0, m=0
      character(len=3), intent(in) :: status


      if (status=="ini") then
         !wir benoetigen hier bereits ausgerechnete hilfsvariablen -> 'getready'

         call CALC_ETOT ! one run for L_totA
         call CALC_ETOT ! another one for L_tot
         L_in_sum = z0; L_out_sum = z0 ! reset boundary terms
         m_in_sum = z0; m_out_sum = z0 ! reset boundary terms

         if (save_tpd /= 0)                    call WRITE_TPD   !Startmodell ins tpd-File

      else if (status=="run") then
         !-----------------------------------------------------------------------
         !  Berechnung von Gesamt-drehimpuls, -masse und -energie
         !-----------------------------------------------------------------------
         ! muss ich immer aufrufen - auch wenn momentan nicht fuer den output benoetigt
         ! andernfalls kann die Summierung an den Raendern nicht stimmen !

         call CALC_ETOT

         if (loglevel >= 2) write(66,"(2(a,ES12.5))") "output.f90:   max. Verdrehung = ",twistmax,"   CFLnumber = ",CFLnumber


         !-----------------------------------------------------------------------
         !  Output waehrend der Rechnung
         !-----------------------------------------------------------------------

         k=k+1; l=l+1; m=m+1

         if (k==save_tpd) then
            call WRITE_TPD
            k = 0
         end if

         if (m==save_dat) then
            call WRITE_DAT
            m = 0
         end if

      else if (status=="fin") then
         ! Ausgabe des letzten Modells

         call WRITE_DAT

      else
         write(66,"(a)")"output.f90:   Aufruf von OUTPUT mit unzulaessigem status-string -> STOP"
         write(* ,"(a)")"output.f90:   Aufruf von OUTPUT mit unzulaessigem status-string -> STOP"
         STOP
      end if

END SUBROUTINE OUTPUT





SUBROUTINE CALC_ETOT
!
! Verwendungszweck:
!  Berechnung Berechnung von Gesamtdrehimpuls, -masse und -energie fuer den Variablensatz X
!
!  Die Vorzeichenkonvention fuer die Randterme ist das Austroemungen aus dem Rechenbereich
!  ein positives Vorzeichen haben. d.h. man muss die Summen der Randterme zu der Summe ueber
!  den Rechenbereich addieren um erhaltene Groessen zu erhalten.
!    also z.B.  m_tot + m_in_sum + m_out_sum = const.
!
!  CALC_ETOT benoetigt einen aktuellen Satz Hilfsvariable (siehe 'use'-liste)
!  Nach einer Konvergenz in NEWTON sind ohnehin alle Hilfsvariable aktuell; am Beginn der
!  Rechnung werden in GETREADY die Hilfsvariable fuer das Startmodell ausgerechnet.
!
!-----------------------------------------------------------------------
! Interne Variable:
!
!   idx       : Schleifenvariable
!   e_int     : Aufsummierung der inneren Energie
!   e_pot     : Aufsummierung der potentiellen Energie
!   e_kin_r   : Aufsummierung der kinetischen Energie mit u_r
!   e_kin_phi : Aufsummierung der kinetischen Energie mit u_phi
!
!   L_in      : Drehimpulstransport ueber den Innenrand waehrend eins Zeitschritts
!   L_out     : Drehimpulstransport ueber den Aussenrand waehrend eins Zeitschritts
!   m_in      : Massentransport ueber den Innenrand waehrend eins Zeitschritts
!   m_out     : Massentransport ueber den Aussenrand waehrend eins Zeitschritts
!
!   twist(np) : Verdrehung durch die Rotation der Scheibe waehrend eines Zeitschrittes, fuer alle Gitterpunkte
!   L_totA    : Gesamtdrehimpuls zum alten Zeitpunkt, zur Berechnung des Drehimpulsfehlers
!
!-----------------------------------------------------------------------

      use config,   only : rkind, np
      use physco ,  only : z0, z1, z8, z12, StBolz, grav, pi
      use primvar,  only : X, XZ, MR, MD, ME, MUr, MUp, MJ
      use global,   only : tst
      use geomvar,  only : S_vol, V_vol, S_flux
      use RBvar,    only : M_star
      use viscvar,  only : Uterm_p, muQ_p
      use advecvar, only : momp_adv, rho_adv
      use eddvar,   only : mu_z

      implicit none

      integer                :: idx
      real(kind=rkind)       :: e_int, e_kin_r, e_kin_phi, e_pot
      real(kind=rkind)       :: L_in, L_out, m_in, m_out
      real(kind=rkind)       :: twist(np)
      real(kind=rkind), save :: L_totA = z1


!     Null-setzen vor der Summierung ueber alle Indizes
      L_tot     = z0
      m_tot     = z0
      e_int     = z0
      e_kin_phi = z0
      lum       = z0
      e_kin_r   = z0
      e_pot     = z0


!-----------------------------------------------------------------------
!    Aufsummierung ueber alle Indizes
!-----------------------------------------------------------------------

      do idx = 3,np-3 ! Summe ueber alle skalaren Zellen

         L_tot     = L_tot     + S_vol(idx) * X(MD,idx) * X(MUp,idx) * z12*( X(MR,idx)  + X(MR,idx+1)  )
         m_tot     = m_tot     + S_vol(idx) * X(MD,idx)
         e_int     = e_int     + S_vol(idx) * X(MD,idx) * X(ME,idx)
         ! Strahlungsenergichte kann ich hier nicht beruecksichtigen da keine z-integrierte Variable

         e_kin_phi = e_kin_phi + z12 * S_vol(idx) * X(MD,idx) * X(MUr,idx)**2

         lum       = lum       + z8*pi * mu_z * X(MJ,idx) * S_vol(idx)

         twist(idx) = X(MUp,idx) * tst / ( pi * ( X(MR,idx)  + X(MR,idx+1) ) )

      end do
      do idx = 3,np-2 ! Summe ueber alle vektoriellen Zellen

         e_kin_r   = e_kin_r   + z12 * V_vol(idx) * X(MD,idx) * X(MUr,idx)**2
         e_pot     = e_pot     - grav * M_star / X(MR,idx) * V_vol(idx) * X(MD,idx)

      end do


!-----------------------------------------------------------------------
!    Randterme
!-----------------------------------------------------------------------

      idx = 3    ! Randterme fuer die innereste skalare Zelle

      L_in      = - S_flux(idx) * XZ(MR,idx) * momp_adv(idx)                                                         & ! advection
                  + z12*pi * muQ_p(idx)   * XZ(MR,idx)**2   * Uterm_p(idx)   * tst                                     ! viscosity
      m_in      = - S_flux(idx)   * rho_adv(idx)                                                                       ! advection

      idx = np-3 ! Randterme fuer die auesserste skalare Zelle

      L_out     = + S_flux(idx+1) * XZ(MR,idx+1) * momp_adv(idx+1)                                                   & ! advection
                  - z12*pi * muQ_p(idx+1) * XZ(MR,idx+1)**2 * Uterm_p(idx+1) * tst                                     ! viscosity
      m_out     = + S_flux(idx+1) * rho_adv(idx+1)                                                                     ! advection

!-----------------------------------------------------------------------
!    Berechnung der integraler Parameter
!-----------------------------------------------------------------------

      e_tot     = e_int + e_kin_r + e_kin_phi + e_pot

      L_err     = abs( (L_tot + L_in + L_out - L_totA) / L_totA )
      L_totA    = L_tot

      T_eff     = sqrt( sqrt( lum / ( sum( S_vol(:) ) * StBolz ) ) )

      ! fuer alle 'full-size' vektoriellen Zellen: idx = 4 .. np-3
      CFLnumber = maxval( abs( S_flux(4:np-3) ) / V_vol(4:np-3) )

      ! fuer alle skalare Zellen: idx = 3 .. np-3
      twistmax  = maxval( abs( twist(3:np-3) ) )

      L_out_sum = L_out_sum + L_out
      L_in_sum  = L_in_sum  + L_in
      m_out_sum = m_out_sum + m_out
      m_in_sum  = m_in_sum  + m_in

      accr_in   =   m_in / tst
      accr_out  =   m_out / tst


END SUBROUTINE CALC_ETOT





SUBROUTINE WRITE_TPD
!
! Verwendungszweck:
!     Schreiben des .tpd Files
!
!-----------------------------------------------------------------------
! Interne Variable:
!
!   datfilelen   : Die Laenge von "datfile"
!   modnrchar    : Die Modellnummer "modnr" als CHARACTER -Variable
!   modnrcharlen : Die Laenge von "modchar"
!   filename     : Name der Output-Files
!   filenamelen  : genutzte Wortlaenge in "filename"
!   archive_reference : Name des zum Kompilat gehoerigen Archiv zip-files
!   backup_reference  : Name des letzten vollstaendigen Backups
!
!-----------------------------------------------------------------------

      use global, only : t, modnr, relax, datfile, logfilename

      implicit none

      integer              :: datfilelen, modnrcharlen, filenamelen
      character(len=9)     :: modnrchar
      character(len=31)    :: archive_reference
      character(len=29)    :: backup_reference
      character(len=112)   :: filename   ! 2 mal 50 +".???????"+".dat"
      logical, save        :: init = .true.

!-----------------------------------------------------------------------
!    Zusammensetzen des File-Namens und oeffnen des tpd-Files
!-----------------------------------------------------------------------
      if (init) then
         datfilelen = len_trim ( datfile )

         if (relax) then
            filename = datfile(1:datfilelen) // ".rel.tpd"
            filenamelen = datfilelen + 8
         else
            write(modnrchar,"(I9)") modnr
            modnrchar = adjustl( modnrchar )
            modnrcharlen = len_trim( modnrchar )

            filename = datfile(1:datfilelen) // "." // modnrchar(1:modnrcharlen) // ".tpd"
            filenamelen = datfilelen + 1 + modnrcharlen + 4
         end if

         include '../../common-code/archive-reference.txt'
         include '../../common-code/backup-reference.txt'


         open(unit=11,file=filename,action="write",status="new",err=999)

         write(11,"(a,a)") "#    corresponding log file     : ", logfilename
         write(11,"(a,a)") "#    corresponding code archive : ", archive_reference
         write(11,"(a,a)") "#    last full backup           : ", backup_reference
         write(11,"(a)")   "#"
         write(11,"(a)")   "# modnr, t, m_tot, m_in_sum, m_out_sum, accr_in, accr_out, L_tot, L_in_sum, L_out_sum, L_err, e_tot, lum, T_eff"

         init = .false.

      end if

!-----------------------------------------------------------------------
!    Ausgabe
!-----------------------------------------------------------------------

      write(11,"(I9,13(1xES23.15E3))")  modnr, t, m_tot, m_in_sum, m_out_sum, accr_in, accr_out, L_tot, L_in_sum, L_out_sum, L_err, e_tot, lum, T_eff

      RETURN


999   continue

      write(66,"(a,a,a)")"output.f90:   Fehler beim Anlegen des Datenfiles: ",filename(1:filenamelen)," -> STOP"
      write(* ,"(a,a,a)")"output.f90:   Fehler beim Anlegen des Datenfiles: ",filename(1:filenamelen)," -> STOP"
      STOP

END SUBROUTINE WRITE_TPD





SUBROUTINE WRITE_DAT
!
! Verwendungszweck:
!     Ausgabe der eines Modells, i.e. kompletter Satz primaerer Variable
!     inklusive Header und einigen zusaetzlichen Datenreihen zum Plotten
!
!-----------------------------------------------------------------------
! Interne Variable:
!
!   idx          : Schleifenindex
!   modnrchar    : Die Modellnummer "modnr" als CHARACTER -Variable
!   modnrcharlen : Die Laenge von "modchar"
!   datfilelen   : Die Laenge von "datfile"
!   filename     : Name der Output-Files
!   filenamelen  : genutzte Wortlaenge in "filename"
!   archive_reference : Name des zum Kompilat gehoerigen Archiv zip-files
!   backup_reference  : Name des letzten vollstaendigen Backups
!
!-----------------------------------------------------------------------

      use config,  only : np
      use global,  only : t, modnr, loglevel, datfile, logfilename, relax
      use primvar, only : nv, X, MR, MD, ME, MUr, MUp, MJ, MH
      use zvar,    only : pseudo2D, H_p, rho0
      use RBvar,   only : M_star, L_star, R_in, rho0_in, e_in, R_out, rho0_out, e_out
      use matvar,  only : Pgas0, Tgas, OPAros, OPAplk
      use physco,  only : z0

      implicit none

      integer                       :: idx, datfilelen, modnrcharlen, filenamelen
      character(len=9)              :: modnrchar
      character(len=31)             :: archive_reference
      character(len=29)             :: backup_reference
      character(len=112)            :: filename   ! 2 mal 50 +".???????"+".dat"

!-----------------------------------------------------------------------
!    Zusammensetzen des File-Namens
!----------------------------------------------------------------------

      datfilelen = len_trim ( datfile )

      if (relax) then
         if (loglevel >= 1) write(66,"(a)") "output.f90:   Ausgabe des des relaxierten Modells"

         filename = datfile(1:datfilelen) // ".0.dat"
         filenamelen = datfilelen + 6

         t = z0 !Modellzeit wird auf 0 zurueckgesetzt

         ! reset boundary values for pseudo cells according to relaxed model
         ! (as the structure of the relaxed model can be vastly different from the initial model)
         rho0_in       = rho0(3)
         X(MD,1:2)     = X(MD,3)

         e_in          = X(ME,3)
         X(ME,1:2)     = X(ME,3)       ! changing the energy upsets the radiation variables, maybe these schould be modified accordingly here as well ?? !debugging

         rho0_out      = rho0(np-3)
         X(MD,np-2:np) = X(MD,np-3)

         e_out         = X(ME,np-3)
         X(ME,np-2:np) = X(ME,np-3)

      else
         if (loglevel >= 1) write(66,"(a,i8,a,ES23.15E3)") "output.f90:   Ausgabe des Modells Nr.",modnr," fuer den Zeitpunkt t = ", t

         write(modnrchar,"(I9)") modnr
         modnrchar = adjustl( modnrchar )
         modnrcharlen = len_trim( modnrchar )

         filename = datfile(1:datfilelen) // "." // modnrchar(1:modnrcharlen) // ".dat"
         filenamelen = datfilelen + 1 + modnrcharlen + 4
      end if

      include '../../common-code/archive-reference.txt'
      include '../../common-code/backup-reference.txt'


!-----------------------------------------------------------------------
!    Ausgabe
!-----------------------------------------------------------------------

      open(unit=10,file=filename,action="write",status="new",err=999)

         write(10,"(a)")                             "# 1DD model file"
         write(10,"(a,a)")                           "#    corresponding log file     : ", logfilename
         write(10,"(a,a)")                           "#    corresponding code archive : ", archive_reference
         write(10,"(a,a)")                           "#    last full backup           : ", backup_reference
         write(10,"(a)")                             "#"
         write(10,"(a)")                             "# np, nv, pseudo2D, t, M_star, L_star, R_in, rho0_in, e_in, R_out, rho0_out, e_out"
         write(10,"(a,2(1xi8),1xL1,9(1xES23.15E3))") "#",np, nv, pseudo2D, t, M_star, L_star, R_in, rho0_in, e_in, R_out, rho0_out, e_out
         write(10,"(a)")                             "#"
         write(10,"(a)")                             "# data format:   1   2   3   4   5    6    7   8  | 9    10     11    12      13"
         write(10,"(a)")                             "#               idx, MR, MD, ME, MUr, MUp, MJ, MH | H_p, Pgas0, Tgas, OPAros, OPAplk"
         do idx = 1,np
            write(10,"(i8,12(1xES23.15E3))") &
            idx, X(MR,idx), X(MD,idx) , X(ME,idx), X(MUr,idx), X(MUp,idx), X(MJ,idx) , X(MH,idx),&
            H_p(idx), Pgas0(idx) , Tgas(idx) , OPAros(idx) , OPAplk(idx)
         end do

      close(unit=10)


      RETURN

999   continue

      write(66,"(a,a,a)")"output.f90:   Fehler beim Anlegen des Datenfiles: ",filename(1:filenamelen)," -> output uebersprungen ! "

END SUBROUTINE WRITE_DAT


END MODULE OUTVAR
